In [ ]:
import geoengine as ge
from geoengine.workflow_builder.operators import GdalSource, TemporalRasterAggregation, RasterStacker, RenameBands, \
Expression,RasterTypeConversion, RasterBandDescriptor
import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from typing import Tuple, Optional
In [ ]:
ge.initialize("https://zentrale.app.geoengine.io/api", token="51c46772-6c55-4867-8e42-c9fd8b614820")
In [ ]:
def _query_rectangle(*,
center: Tuple[float, float],
time: np.datetime64,
radius_px: float = 512) -> ge.QueryRectangle:
resolution = ge.SpatialResolution(10, 10)
bbox = ge.BoundingBox2D(
xmin=center[0] - resolution.x_resolution * radius_px,
xmax=center[0] + resolution.x_resolution * radius_px,
ymin=center[1] - resolution.y_resolution * radius_px,
ymax=center[1] + resolution.y_resolution * radius_px,
)
return ge.QueryRectangle(
spatial_bounds=bbox,
time_interval=ge.TimeInterval(time),
srs='EPSG:32632',
resolution=resolution,
)
async def compute_rgb_array(*,
center: Tuple[float, float],
time: np.datetime64,
radius_px: float = 512,
) -> xr.DataArray:
'''Compute RGB of area'''
def cloud_free(band: GdalSource):
return Expression(
expression = 'if B > 2 && B < 8 { A } else { NODATA }',
source=RasterStacker(
sources=[
band,
RasterTypeConversion(
GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:SCL`"),
output_data_type='U16',
),
],
)
)
query_rectangle = _query_rectangle(
center=center,
time=time,
radius_px=radius_px,
)
workflow = TemporalRasterAggregation(
aggregation_type='mean', # 'percentileEstimate',
# percentile=0.5,
granularity='months',
window_size=1,
ignore_no_data=True,
source=RasterStacker(
sources=[
cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B04`")),
cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B03`")),
cloud_free(GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B02`")),
],
rename=RenameBands.rename(['red', 'green', 'blue']),
)
)
workflow = ge.register_workflow(workflow)
data_array = await workflow.raster_stream_into_xarray(
query_rectangle=query_rectangle,
clip_to_query_rectangle=True,
bands=[0,1,2], # TODO: improve for user, default = all? where are the band names?
)
return data_array
async def compute_water_bodies(*,
center: Tuple[float, float],
time: np.datetime64,
radius_px: float = 512,
water_body_time_aggregate_months: int = 1,
) -> xr.DataArray:
'''Compute water body classification'''
query_rectangle = _query_rectangle(
center=center,
time=time,
radius_px=radius_px,
)
workflow = TemporalRasterAggregation(
aggregation_type='max',
granularity='months',
window_size=water_body_time_aggregate_months,
ignore_no_data=True,
source=Expression(
expression='''
let ndwi = (A-B)/(A+B);
if C <= 2 || C >= 8 {
NODATA
} else if ndwi >= 0.2 {
3
} else if ndwi >= 0 {
2
} else if ndwi >= -0.3 {
1
} else {
0
}
''',
output_type='U8',
# output_band=RasterBandDescriptor( # TODO: this is broken in the backend
# "water bodies",
# measurement=ge.ClassificationMeasurement(
# "water bodies",
# classes={
# 0: 'Drought, non-aqueous surfaces',
# 1: 'Moderate drought, non-aqueous surfaces',
# 2: 'Flooding, humidity',
# 3: 'Water surface',
# }
# ),
# ),
source=
RasterStacker(
sources=[
GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B03`"),
GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:B08`"),
RasterTypeConversion(
GdalSource("_:5779494c-f3a2-48b3-8a2d-5fbba8c5b6c5:`UTM32N:SCL`"),
output_data_type='U16',
),
],
),
),
)
workflow = ge.register_workflow(workflow)
data_array = await workflow.raster_stream_into_xarray(
query_rectangle=query_rectangle,
clip_to_query_rectangle=True,
bands=[0], # TODO: improve for user, default = all? where are the band names?
)
return data_array
def plot(rbg_array: xr.DataArray, water_bodies_array: Optional[xr.DataArray] = None, height_px: int = 1024):
px = 1/plt.rcParams['figure.dpi'] # pixel in inches
fig, axes = plt.subplot_mosaic("AB")
fig.set_figheight(height_px*px)
fig.set_figwidth(2 * height_px*px)
rbg_array.isel(time=0).plot.imshow(
rgb="band",
vmax=4000,
ax=axes['A'],
)
if water_bodies_array is None:
fig.show()
return
water_bodies_array.isel(band=0, time=0, drop=True).plot.imshow(
colors=[plt.cm.Paired.colors[i] for i in [2, 3, 0, 1]],
levels=5,
vmin=0,
vmax=4,
ax=axes['B']
)
axes['B'].images[-1].colorbar.set_ticks(
[0,1,2,3],
labels=[
'Drought, non-aqueous surfaces',
'Moderate drought, non-aqueous surfaces',
'Flooding, humidity',
'Water surface',
],
)
fig.show()
fig, axes = plt.subplot_mosaic("C")
fig.set_figheight(height_px*px)
fig.set_figwidth(height_px*px)
rbg_array.isel(time=0).plot.imshow(
rgb="band",
vmax=4000,
ax=axes['C'],
)
water_bodies_array.isel(band=0, time=0, drop=True).plot.imshow(
colors=[(0,0,0,0), plt.cm.Paired.colors[1]],
levels=2,
vmin=3,
vmax=4,
add_colorbar=False,
ax=axes['C']
)
fig.show()
In [ ]:
marburg_center_utm = [483843, 5627614]
In [ ]:
marburg_rgb = await compute_rgb_array(center=marburg_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn(
In [ ]:
marburg_water = await compute_water_bodies(center=marburg_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn(
In [ ]:
plot(rbg_array=marburg_rgb, water_bodies_array=marburg_water)
/tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show() /tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show() /home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast xx = (xx * 255).astype(np.uint8)
In [ ]:
weimar_center_utm = [479763, 5623423]
weimar_rgb = await compute_rgb_array(center=weimar_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
weimar_water = await compute_water_bodies(center=weimar_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
plot(rbg_array=weimar_rgb, water_bodies_array=weimar_water)
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn( /home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn( /tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show() /tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show() /home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/matplotlib/cm.py:478: RuntimeWarning: invalid value encountered in cast xx = (xx * 255).astype(np.uint8)
In [ ]:
koeln_center_utm = [356766, 5644819]
koeln_rgb = await compute_rgb_array(center=koeln_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
koeln_water = await compute_water_bodies(center=koeln_center_utm, time=np.datetime64("2022-07-01T00:00:00"))
plot(rbg_array=koeln_rgb, water_bodies_array=koeln_water)
/home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn( /home/beilschmidt/git/geoengine-python/env/lib/python3.10/site-packages/rasterio/windows.py:314: RasterioDeprecationWarning: The height, width, and precision parameters are unused, deprecated, and will be removed in 2.0.0. warnings.warn( /tmp/ipykernel_9616/3560293408.py:180: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show() /tmp/ipykernel_9616/3560293408.py:202: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure. fig.show()
